

use "stata data/A1 Standard regression dataset", clear
capture drop _merge
keep if inrange(bin,-8,8)

keep if inrange(hour,14,17)

keep if event_day == 1

capture drop _merge

rename temp_f_mean temp_f
gen temp_f_2 = temp_f*temp_f
gen ln_kwh = ln(kwh)

label variable temp_f "Temperature (F)"
label variable temp_f_2 "Temperature (F) squared"
label variable pdp "2-6pm reductions"
label variable eligible "Eligible for PDP"
gen dow = dow(date)


gen month = month(date)


drop if total > 50000
drop if total < 4000


gen coastal = 0
replace coastal = 1 if inlist(climate,5,6,2)


drop if kwh == 0

save "stata data/main regs dataset with optionally", replace

drop if optionally_enrolled == 1

preserve
use "temp/opt out new extract updated", clear
drop if action_status == "Closed SA" //do this so I"m not counting the people that just seem to drop out

gen dropout_date_ym = mofd(pdp_end)
gen dropout_during_first_summer = 0
replace dropout_during_first_summer  = 1 if ///
	inrange(dropout_date_ym,665,669) & dropout_date_ym !=.
keep sa pdp_end
save "stata data/drop out test", replace
restore

merge m:1 sa using "stata data/drop out test"
drop if _merge == 2
order pdp_end

//this makes sure I've got things turning on and off at the right time. very few changes
replace pdp = 0 if date >pdp_end & year == 2015 & pdp_end !=. & pdp_end <mdy(11,1,2015)
replace pdp = 1 if date <= pdp_end & year == 2015 & pdp_end !=. & pdp_end <mdy(11,1,2015)


replace temp_f = . if temp_f > 130

egen hour_of_sample = group(hour date)

save "stata data/main regs dataset", replace

//small diversion to get the list of sas that are in regressions
use "stata data/main regs dataset", clear
keep sa
gduplicates drop
save "stata data/sas_used_in_regressions", replace

use "stata data/main regs dataset", clear



///Table 2 regressions -->
//first stage

ivreghdfe pdp eligible temp_f temp_f_2  if event_day == 1 , absorb(sa#hour#dow  date#hour) cluster( sa hour_of_sample) 
unique sa 
estadd scalar firms = `r(unique)'

sum temp_f if event_day == 1
estadd scalar temp_f= `r(mean)'


regsave using "stata data/regsave/A1_new_first_stage", replace ci t p ///
addlabel(type,"event",cz,"all",firms,`e(firms)', avg_kwh,`e(avg_kwh)', f,`e(F)', obs, `e(N)', temp_f, `e(temp_f)')


ivreghdfe pdp eligible temp_f temp_f_2  if event_day == 1 & coastal == 1 , absorb(sa#hour#dow  date#hour) cluster( sa hour_of_sample) 
unique sa if coastal == 1
estadd scalar firms = `r(unique)'

sum temp_f if coastal == 1
estadd scalar temp_f= `r(mean)'

regsave using "stata data/regsave/A1_new_coastal_first_stage", replace ci t p ///
addlabel(type,"event",cz,"coastal",firms,`e(firms)', avg_kwh,`e(avg_kwh)', f,`e(F)' , obs, `e(N)', temp_f, `e(temp_f)')


ivreghdfe pdp eligible temp_f temp_f_2  if event_day == 1 & coastal == 0 , absorb(sa#hour#dow  date#hour) cluster( sa hour_of_sample) 


unique sa if coastal == 0
estadd scalar firms = `r(unique)'

sum temp_f if coastal == 0
estadd scalar temp_f= `r(mean)'


regsave using "stata data/regsave/A1_new_inland_first_stage", replace ci t p ///
addlabel(type,"event",cz,"inland",firms,`e(firms)', avg_kwh,`e(avg_kwh)', f,`e(F)' , obs, `e(N)', temp_f, `e(temp_f)')

///<----Table 2 regressions 


///----> Table 3 output 




ivreghdfe ln_kwh temp_f temp_f_2 (pdp=eligible)  if event_day == 1 , absorb(sa#hour#dow  date#hour) cluster(sa hour_of_sample) 

unique sa 
estadd scalar firms = `r(unique)'


sum kwh if event_day == 1 & eligible == 0 
estadd scalar avg_kwh = `r(mean)'

sum temp_f if event_day == 1 
estadd scalar temp_f= `r(mean)'


regsave using "stata data/regsave/FE_new_all_event", replace ci t p ///
addlabel(type,"event",cz,"all",firms,`e(firms)', avg_kwh,`e(avg_kwh)', temp_f, `e(temp_f)', clusters,`e(N_clust)', obs, `e(N)' )


ivreghdfe ln_kwh temp_f temp_f_2 (pdp=eligible)  if event_day == 1 & coastal == 1  , absorb(sa#hour#dow  date#hour) cluster( sa hour_of_sample) 
unique sa if coastal == 1
estadd scalar firms = `r(unique)'


sum kwh if event_day == 1 & coastal == 1 & eligible == 0
estadd scalar avg_kwh = `r(mean)' 

sum temp_f if event_day == 1 & coastal == 1
estadd scalar temp_f= `r(mean)'

regsave using "stata data/regsave/FE_new_all_coastal_event", replace ci t p ///
addlabel(type,"event",cz,"coastal",firms,`e(firms)', avg_kwh,`e(avg_kwh)', temp_f, `e(temp_f)', clusters,`e(N_clust)', obs, `e(N)' )


ivreghdfe ln_kwh temp_f temp_f_2 (pdp=eligible)  if event_day == 1 & coastal == 0  , absorb(sa#hour#dow  date#hour) cluster( sa hour_of_sample) 
unique sa if coastal == 0 
estadd scalar firms = `r(unique)'


sum kwh if event_day == 1 & coastal == 0 & eligible == 0
estadd scalar avg_kwh = `r(mean)' 

sum temp_f if event_day == 1 & coastal == 0
estadd scalar temp_f= `r(mean)'

regsave using "stata data/regsave/FE_new_all_inland_event", replace ci t p ///
addlabel(type,"event",cz,"inland",firms,`e(firms)', avg_kwh,`e(avg_kwh)', temp_f, `e(temp_f)', clusters,`e(N_clust)', obs, `e(N)' )





//// Table 4 regs--->
///
///Non-event days
///

use if inrange(bin,-8,8) using "stata data/A1 Standard regression dataset", clear

keep if event_day == 0

replace temp_f = . if temp_f > 130

capture drop _merge

rename temp_f_mean temp_f
gen temp_f_2 = temp_f*temp_f
gen ln_kwh = ln(kwh)

gen dow = dow(date)
gen month = month(date)

drop if total > 50000
drop if total < 4000
gen coastal = 0
replace coastal = 1 if inlist(climate,5,6,2)

drop if kwh == 0

drop if optionally_enrolled == 1

egen hour_of_sample = group(hour date)

rename stn_call stn_call_str
encode stn_call_str, gen(stn_call)


ivreghdfe ln_kwh temp_f temp_f_2 (pdp=eligible)  if event_day == 0 , absorb(sa#hour#dow  date#hour) cluster(sa hour_of_sample)
unique sa 
estadd scalar firms = `r(unique)'


sum kwh if event_day == 0 & eligible == 0 
estadd scalar avg_kwh = `r(mean)'

sum temp_f if event_day == 0 
estadd scalar temp_f= `r(mean)'


regsave using "stata data/regsave/FE_new_all_non_event", replace ci t p ///
addlabel(type,"event",cz,"all",firms,`e(firms)', avg_kwh,`e(avg_kwh)', temp_f, `e(temp_f)', clusters,`e(N_clust)', obs, `e(N)' )


ivreghdfe ln_kwh temp_f temp_f_2 (pdp=eligible)  if event_day == 0 & coastal == 1  , absorb(sa#hour#dow  date#hour) cluster( sa hour_of_sample) 
unique sa if coastal == 1
estadd scalar firms = `r(unique)'


sum kwh if event_day == 0 & coastal == 1 & eligible == 0
estadd scalar avg_kwh = `r(mean)' 

sum temp_f if event_day == 0 & coastal == 1
estadd scalar temp_f= `r(mean)'

regsave using "stata data/regsave/FE_new_all_coastal_non_event", replace ci t p ///
addlabel(type,"event",cz,"coastal",firms,`e(firms)', avg_kwh,`e(avg_kwh)', temp_f, `e(temp_f)', clusters,`e(N_clust)', obs, `e(N)' )

ivreghdfe ln_kwh temp_f temp_f_2 (pdp=eligible)  if event_day == 0 & coastal == 0  , absorb(sa#hour#dow  date#hour) cluster( sa hour_of_sample) 
unique sa if coastal == 0 
estadd scalar firms = `r(unique)'


sum kwh if event_day == 0 & coastal == 0 & eligible == 0
estadd scalar avg_kwh = `r(mean)' 

sum temp_f if event_day == 0 & coastal == 0
estadd scalar temp_f= `r(mean)'

regsave using "stata data/regsave/FE_new_all_inland_non_event", replace ci t p ///
addlabel(type,"event",cz,"inland",firms,`e(firms)', avg_kwh,`e(avg_kwh)', temp_f, `e(temp_f)', clusters,`e(N_clust)', obs, `e(N)' )

///<--- end table 4 regs




/// Table 5 regs -->

///
///Bins based on individual temps (5/24/21) --->
///
use "stata data/main regs dataset", clear

capture egen hour_of_sample = group(hour date)

capture drop temp_bin
egen temp_bin = cut(temp_f) ,at(70,80,90,100)
replace temp_bin = 100 if temp_f >= 100
replace temp_bin = 70 if temp_f <= 70

tab temp_bin if year ==2015 & coastal == 0
tab temp_bin if year ==2015 & coastal == 0 & pdp == 1
gunique sa if temp_bin==100 & year ==2015 & coastal == 0
gunique sa if temp_bin==100 & year ==2015 & coastal == 0 & eligible == 1
gunique sa if temp_bin==100 & year ==2015 & coastal == 0 & pdp == 1


///just if I want a count
preserve
capture drop counter
gen counter = 1
collapse (sum) counter, by(coastal temp_bin)
save "stata data/obs_counts_by_temp_bin", replace
restore

tab temp_bin, gen(temp_bin_fe)


foreach num of numlist 1/4 {
gen pdp_by_temp_bin_`num' = pdp * temp_bin_fe`num'
gen eligible_by_temp_bin_`num' = eligible * temp_bin_fe`num'
}

ivreghdfe ln_kwh temp_f temp_f_2 (pdp_by_temp_bin_*  = eligible_by_temp_bin_* )  if event_day == 1  , absorb(sa#hour#dow  date#hour) cluster( sa hour_of_sample) 

unique sa 
estadd scalar firms = `r(unique)'


sum kwh if event_day == 1 & eligible == 0 
estadd scalar avg_kwh = `r(mean)'

sum temp_f if event_day == 1 
estadd scalar temp_f= `r(mean)'


regsave using "stata data/regsave/FE_new_all_individual_tempbins_event", replace ci t p ///
addlabel(type,"event",cz,"all",firms,`e(firms)', avg_kwh,`e(avg_kwh)', temp_f, `e(temp_f)', clusters,`e(N_clust)', obs, `e(N)' )


ivreghdfe ln_kwh temp_f temp_f_2 (pdp_by_temp_bin_*  = eligible_by_temp_bin_* )   if event_day == 1 & coastal == 1  , absorb(sa#hour#dow  date#hour) cluster( sa hour_of_sample) 
unique sa if coastal == 1
estadd scalar firms = `r(unique)'

sum kwh if event_day == 1 & coastal == 1 & eligible == 0
estadd scalar avg_kwh = `r(mean)' 

sum temp_f if event_day == 1 & coastal == 1
estadd scalar temp_f= `r(mean)'

regsave using "stata data/regsave/FE_new_all_individual_tempbins_coastal_event", replace ci t p ///
addlabel(type,"event",cz,"coastal",firms,`e(firms)', avg_kwh,`e(avg_kwh)', temp_f, `e(temp_f)', clusters,`e(N_clust)', obs, `e(N)' )



ivreghdfe ln_kwh temp_f temp_f_2 (pdp_by_temp_bin_*  = eligible_by_temp_bin_* )   if event_day == 1 & coastal == 0  , absorb(sa#hour#dow  date#hour) cluster( sa hour_of_sample) 
unique sa if coastal == 0 
estadd scalar firms = `r(unique)'


sum kwh if event_day == 1 & coastal == 0 & eligible == 0
estadd scalar avg_kwh = `r(mean)' 

sum temp_f if event_day == 1 & coastal == 0
estadd scalar temp_f= `r(mean)'

regsave using "stata data/regsave/FE_new_all_individual_tempbins_inland_event", replace ci t p ///
addlabel(type,"event",cz,"inland",firms,`e(firms)', avg_kwh,`e(avg_kwh)', temp_f, `e(temp_f)', clusters,`e(N_clust)', obs, `e(N)' )


///
///Bins based on individual temps <---
///
/// Table 5 regs <--



/// Table 6 regs -->

///
///Bins based on forecasted temps
///
use "stata data/main regs dataset", clear

drop _merge
merge m:1 date using "raw data/trigger and forcasted temps" //this file was provided to me by PG&E 
keep if _merge == 3


gen temp_bin = .
replace temp_bin = 97 if inrange(forecast_temps,96,97)
replace temp_bin = 99 if inrange(forecast_temps,98,99)
replace temp_bin = 101 if inrange(forecast_temps,100,101)
replace temp_bin = 102 if inrange(forecast_temps,102,104)

tab temp_bin, gen(temp_bin_fe)

foreach num of numlist 1/4 {
gen pdp_by_temp_bin_`num' = pdp * temp_bin_fe`num'
gen eligible_by_temp_bin_`num' = eligible * temp_bin_fe`num'
}




ivreghdfe ln_kwh temp_f temp_f_2 (pdp_by_temp_bin_*  = eligible_by_temp_bin_* )  if event_day == 1  , absorb(sa#hour#dow  date#hour) cluster( sa hour_of_sample) 


unique sa 
estadd scalar firms = `r(unique)'


sum kwh if event_day == 1 & eligible == 0 
estadd scalar avg_kwh = `r(mean)'

sum temp_f if event_day == 1 
estadd scalar temp_f= `r(mean)'


regsave using "stata data/regsave/FE_new_all_tempbins_event", replace ci t p ///
addlabel(type,"event",cz,"all",firms,`e(firms)', avg_kwh,`e(avg_kwh)', temp_f, `e(temp_f)', clusters,`e(N_clust)', obs, `e(N)' )



ivreghdfe ln_kwh temp_f temp_f_2 (pdp_by_temp_bin_*  = eligible_by_temp_bin_* )   if event_day == 1 & coastal == 1  , absorb(sa#hour#dow  date#hour) cluster( sa hour_of_sample) 
unique sa if coastal == 1
estadd scalar firms = `r(unique)'


sum kwh if event_day == 1 & coastal == 1 & eligible == 0
estadd scalar avg_kwh = `r(mean)' 

sum temp_f if event_day == 1 & coastal == 1
estadd scalar temp_f= `r(mean)'

regsave using "stata data/regsave/FE_new_all_tempbins_coastal_event", replace ci t p ///
addlabel(type,"event",cz,"coastal",firms,`e(firms)', avg_kwh,`e(avg_kwh)', temp_f, `e(temp_f)', clusters,`e(N_clust)', obs, `e(N)' )




ivreghdfe ln_kwh temp_f temp_f_2 (pdp_by_temp_bin_*  = eligible_by_temp_bin_* )   if event_day == 1 & coastal == 0  , absorb(sa#hour#dow  date#hour) cluster( sa hour_of_sample) 
unique sa if coastal == 0 
estadd scalar firms = `r(unique)'


sum kwh if event_day == 1 & coastal == 0 & eligible == 0
estadd scalar avg_kwh = `r(mean)' 

sum temp_f if event_day == 1 & coastal == 0
estadd scalar temp_f= `r(mean)'

regsave using "stata data/regsave/FE_new_all_tempbins_inland_event", replace ci t p ///
addlabel(type,"event",cz,"inland",firms,`e(firms)', avg_kwh,`e(avg_kwh)', temp_f, `e(temp_f)', clusters,`e(N_clust)', obs, `e(N)' )




///
///Bins based on forecasted temps <---
///
/// Table 6 regs <--



//// Table 7 regs --> 

///
///Naics codes regressions
///

*customer facing from
*http://www.bls.gov/iag/tgs/iag07.htm


use "stata data/main regs dataset", clear


rename stn_call stn_call_str
encode stn_call_str, gen(stn_call)


capture drop test
gen customer_facing = 0
replace customer_facing = 1 if inlist(naics_most_2,44,45,71,72, 62)

gen no_customer_facing = 0
replace no_customer_facing = 1 if customer_facing == 0


keep if naics_most_2 !=.

foreach customer of numlist 0/1 {


preserve
keep if customer == `customer'

ivreghdfe ln_kwh temp_f temp_f_2 (pdp=eligible)  if event_day == 1 , absorb(sa#hour#dow  date#hour) cluster( sa hour_of_sample) 
unique sa 
estadd scalar firms = `r(unique)'


sum kwh if event_day == 1 
estadd scalar avg_kwh = `r(mean)'

sum temp_f if event_day == 1 
estadd scalar temp_f= `r(mean)'


regsave using "stata data/regsave/FE_new_all_customer_`customer'", replace ci t p ///
addlabel(type,"event",cz,"all",firms,`e(firms)', avg_kwh,`e(avg_kwh)', temp_f, `e(temp_f)', clusters,`e(N_clust)', obs, `e(N)',customer_facing,`customer' )




ivreghdfe ln_kwh temp_f temp_f_2 (pdp=eligible)  if event_day == 1 & coastal == 1 , absorb(sa#hour#dow  date#hour) cluster( sa hour_of_sample) 
unique sa if coastal == 1
estadd scalar firms = `r(unique)'


sum kwh if event_day == 1 & coastal == 1
estadd scalar avg_kwh = `r(mean)' 

sum temp_f if event_day == 1 & coastal == 1
estadd scalar temp_f= `r(mean)'

regsave using "stata data/regsave/FE_new_all_coastal_customer_`customer'", replace ci t p ///
addlabel(type,"event",cz,"coastal",firms,`e(firms)', avg_kwh,`e(avg_kwh)', temp_f, `e(temp_f)', clusters,`e(N_clust)', obs, `e(N)' ,customer_facing,`customer' )




ivreghdfe ln_kwh temp_f temp_f_2 (pdp=eligible)  if event_day == 1 & coastal == 0  , absorb(sa#hour#dow  date#hour) cluster( sa hour_of_sample) 
unique sa if coastal == 0
estadd scalar firms = `r(unique)'


sum kwh if event_day == 1 & coastal == 0
estadd scalar avg_kwh = `r(mean)' 

sum temp_f if event_day == 1 & coastal == 0
estadd scalar temp_f= `r(mean)'

regsave using "stata data/regsave/FE_new_all_inland_customer_`customer'", replace ci t p ///
addlabel(type,"event",cz,"inland",firms,`e(firms)', avg_kwh,`e(avg_kwh)', temp_f, `e(temp_f)', clusters,`e(N_clust)', obs, `e(N)' ,customer_facing,`customer' )

restore
}
*
///<---end table 7 regs




///---> Figure 5 regression

use "stata data/A1 data -8 to 8 HUGE", clear

merge m:1 sa using "stata data/sas_used_in_regressions" 
keep if _merge == 3
drop _merge


gen dow = dow(date)
keep if inrange(dow,1,5)


gen eligible = 0
replace eligible = 1 if bin > 0 & year == 2015

compress *



save "stata data/all hours standard regression setup", replace

use "stata data/all hours standard regression setup", clear

replace pdp = 0 if year == 2014


tab hour, gen(hour_fe)

replace pdp = 0 if year == 2014 

foreach num of numlist 1/24{
local  num_m1 = `num' -1
rename hour_fe`num' hour_fe`num_m1'
}


foreach num of numlist 0/23{
gen byte eligible_by_hr_`num' = hour_fe`num'*eligible
gen byte pdp_by_hr_`num' = hour_fe`num'*pdp
}
drop hour_fe*


save "temp/ALL HOURS A1 Standard regression dataset 1-22-2016", replace

use "temp/ALL HOURS A1 Standard regression dataset 1-22-2016", clear

gen ln_kwh = ln(kwh)

keep if event_day == 1

gen temp_f_2 = temp_f^2
gen month = month(date)

gen coastal = 0
replace coastal = 1 if inlist(climate,5,6,2)

keep if coastal == 0

ivreghdfe ln_kwh temp_f temp_f_2  (pdp_by_hr_* =  eligible_by_hr_*) if coastal == 0 , absorb(sa#dow#hour date#hour) cluster(group_dup) 
regsave using "stata data/regsave/A1_ALLHOURS_inland_event", replace ci t p 





*below used for table 9

use "stata data/main regs dataset", clear

capture egen hour_of_sample = group(hour date)

rename stn_call stn_call_str
encode stn_call_str, gen(stn_call)



gen temp_f_old = temp_f
replace temp_f = temp_f - 75
replace temp_f_2 = temp_f^2

capture gen t_pdp = temp_f * pdp
capture gen t_eligible = temp_f* eligible

//inland
ivreghdfe ln_kwh temp_f temp_f_2 (t_pdp pdp = t_eligible eligible)  if event_day == 1 & coastal ==0  , absorb(sa#hour#dow#month  date#hour) cluster( sa hour_of_sample) 
estimates store temp_interaction_inland

unique sa if coastal == 0
estadd scalar firms = `r(unique)'

sum kwh if event_day == 1 & coastal == 0
estadd scalar avg_kwh = `r(mean)' 

sum temp_f_old if event_day == 1 & coastal == 0
estadd scalar temp_f= `r(mean)'

regsave using "stata data/regsave/FE_new_all_inland_temp_f_interact", replace ci t p ///
addlabel(type,"event",cz,"inland",firms,`e(firms)', avg_kwh,`e(avg_kwh)', temp_f, `e(temp_f)', clusters,`e(N_clust)', obs, `e(N)' )





